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ABSTRACT 

In  this  paper  a  relationship  between  two  approaches  towards  construction  of  genuinely 
two-dimensional  upwind  advection  schemes  is  established.  One  of  these  approaches  is  of 
the  control  volume  type  applicable  on  structured  cartesian  meshes.  It  resulted  (see  [14], 
[15])  in  the  compact  high  resolution  schemes  capable  of  maintaining  second  order  accuracy 
in  both  homogeneous  and  inhomogeneous  cases.  Another  one  is  the  fluctuation  splitting 
approach  (see  [11],  [3],  [12],  [17]),  which  is  well  suited  for  triangular  (and  possibly)  un¬ 
structured  meshes.  Understanding  the  relationship  between  these  two  approaches  allows  us 
to  formulate  here  a  new  fluctuation  splitting  high  resolution  (i.e.  possible  use  of  artificial 
compression,  while  maintaining  positivity  property)  scheme.  This  scheme  is  shown  to  be 
linearity  preserving  in  inhomogeneous  as  well  as  homogeneous  cases. 
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1.  Introduction 


The  perfect  scheme  for  numerical  advection  still  awaits  discovery.  Its  attributes  should  include 
high  accuracy  in  regions  of  smooth  flow,  together  with  some  positivity  property  that  avoids  over- 
and  undershoots  near  discontinuities.  It  is  well  known  that  these  two  properties  are  not  both 
attainable  within  the  class  of  linear  schemes.  This  implies  that  successful  schemes  must  be  non¬ 
linear,  (i.e.  data  dependent). 

Thus  the  scheme  must  be  furnished  with  a  monitor  function  that  measures  the  local  smoothness 
of  the  data,  for  example  by  comparing  consecutive  gradients,  as  in  the  flux-limiter  approach.  This 
enlarges  the  stencil  of  the  scheme  beyond  the  minimum  needed  to  attain  a  given  formal  accuracy. 
Such  enlargement  is  undesirable  for  parallel  implementation,  and  seems  to  impair  the  convergence 
of  multigrid  methods,  as  well  as  requiring  supplementary  conditions  close  to  boundaries. 

In  [15]  a  monitor  function  was  introduced  that  compares  gradients  in  different  directions.  This  led 
to  an  advection  scheme  of  the  control-volume  type  with  a  minimally-enlarged  stencil,  and  with  good 
multigrid  capabilities.  Recently,  there  has  also  been  developed  the  class  of  “fluctuation-splitting” 
schemes.  These  schemes  are  intended  to  be  coded  as  a  loop  over  triangular  (tetrahedral)  elements, 
employing  no  data  external  to  the  triangle.  Several  non-linear  variants  of  this  technique  have  been 
devised  that  preserve  positivity. 

In  this  paper  we  show  a  strong  relationship  between  the  two  types  of  schemes.  In  fact,  for  linear 
advection  on  regular  grids,  certain  schemes  in  each  class  turn  out  to  be  identical.  This  enables  the 
transfer  of  insights  and  techniques  from  one  class  to  the  other. 

§2  presents  some  basic  definitions  and  principles  regarding  the  discretization  of  a  conservation 
law.  §3  documents  several  schemes  of  the  control  volume  type,  and  shows  how  limiter  functions 
can  be  introduced  into  them.  §4  briefly  describes  the  fluctuation-splitting  approach,  with  some 
of  the  previously-employed  positivity  devices.  §5  reformulates  the  nonlinear  fluctuation-splitting 
schemes  as  limiter  schemes  and  demonstrates  their  identity  with  some  control-volume  schemes.  The 
performance  of  the  fluctuation-splitting  schemes  is  illustrated  by  some  numerical  experiments  in 
§6.  The  conclusions  of  this  work  are  presented  in  §7. 

The  truly  two-dimensional  control-volume  advection  scheme  that  is  capable  of  producing  second 
order  accurate  steady-state  solutions  for  the  case  of  a  non-zero  source  term  is  formulated  follow¬ 
ing  [15]  in  Appendix  A.  The  high-resolution  fluctuation-splitting  counterpart  of  this  scheme  is 
constructed  in  Appendix  B. 

2.  Conservation  law  and  its  discretization 
Consider  the  scalar  conservation  law 

(2.1)  ut  +  'V-f  =  s, 

where  ^  =  {f,g),  and  its  nonconservative  form,  in  which  A  =  {fu,9u) 

(2.2)  ut  +  X-'Vu  =  s. 

In  what  follows,  unless  otherwise  stated,  we  take  A  =  (a,  b),  with  a,  b  constant 

(2.3)  {dt  +  adx  +  bdy)u  =  s. 

This  assumption  does  not  impair  the  generality  of  the  construction.  This  is  because  using  the 
conservative  linearization  procedure  developed  in  [2]  general  nonlinear  conservation  laws  can  be 
represented  locally  by  a  linear  constant  coefficient  advection  equation  of  the  form  (2.3). 

We  are  interested  here  in  the  steady-state  solutions  (2.1),  i.e.  the  case  when  dtu  —  0. 

It  will  appear  very  useful  to  note  that  a  steady-state  solution  of  Eq.(2.1)  in  the  homogeneous  case 
(s  =  0,  the  case  considered  through  the  entire  paper,  except  for  the  Appendix)  has  the  following 
property:  it  is  constant  along  the  direction  of  the  convection  velocity  (characteristic  direction). 
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Figure  1.  Control  volume. 


The  time  stepping  (forward  Euler)  procedure  is  considered  here  to  be  only  a  means  to  reach  the 
steady-state.  The  solution  update  formula  can  be  given  by  the  following 

(2.4) 

k 

i.e.  the  solution  value  in  each  grid  point  at  the  new  time  level  can  be  represented  as  a  combination 
of  the  solution  values  from  the  previous  time  level. 

Applying  the  TVD  concept  to  characterize  the  schemes  producing  non-oscillatory  solutions  ap¬ 
pears  to  be  too  restrictive  in  two  dimensions  (see  [4]).  Therefore,  we  shall  follow  Spekreijse  [16] 
and  use  the  concept  of  positivity. 

Definition  2.1.  A  scheme  is  said  to  be  of  the  positive  type  if  the  solution  update  formula  can  be 
written  in  the  form  (2.4)  in  such  a  way  that  Ck  >  0,  Vfc. 

Solutions  obtained  by  the  means  of  positive  schemes  will  satisfy  a  certain  maximum  principle. 


3.  Control  volume  approach 

Assume  that  the  computational  domain  is  divided  into  square  cells  (control  volumes)  each  one 
associated  with  the  cell-averaged  value  of  the  solution  approximation  (see  Fig.l).  Integrating 
Eq.(2.1)  over  a  control  volume  and  applying  Gauss  theorem  we  obtain 


/  ut  dxdy  -1-  ® 

I  Co  Jdc 


uX  ■  dn 


-IL 


s  dxdy 


Since  our  purpose  is  to  construct  a  second  order  accurate  advection  scheme,  it  is  suIRcient  to 
approximate  Eq.(3.1)  using  one-  and  two-  dimensional  mid-point  quadrature  rules.  Thus,  we  obtain 


=  '^0  -  -J-U+  -  f-  +  9+  -  g-]  +  So  At, 


where  /+,  f-,g+,  g-  are  numerical  fluxes  computed  at  the  centers  of  each  side  of  the  control  volume 
Cq.  The  question  is  how  to  compute  them  in  order  to  obtain  an  advection  scheme  with  desired 
properties. 

In  order  to  simplify  our  further  discussion  in  this  section  we  assume  that 


0  <  a  <  6 
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Figure  2.  Stencils  of  the  first  order  schemes:  a)  dimensional  upwinding:  b)  N  scheme. 


3.1.  Dimensional  upwinding.  The  simplest  example  of  a  positive  linear  scheme  is  given  by  the 
dimensionally  upwinded  scheme  U  with  the  following  fluxes 


(3.4) 


=  auz 
9-  =  bus 


This  corresponds  to  the  following  update  formula 

(3.5)  =Uq  -  ^[b{uo  -  Us)  +  a{uo  -  tis)] 

(Fig.2(a)  presents  the  corresponding  stencil).  However,  this  scheme  suffers  from  an  excessive  nu¬ 
merical  diffusion. 


3.2.  Zero  cross-diffusion  schemes.  Here  we  present  the  zero  cross-diffusion  2D  scheme  that  was 
used  in  [14], [15]  for  constructing  a  nonlinear  zero  cross  diffusion  positive  scheme  (see  also  §3.4  and 
§3.5).  The  fluxes  are 

=  au3  -  lb{u3  - 
=  bus  -  \a{us  - 

These  flux  formulae  can  be  motivated  by  considering  the  characteristics  dy/dx  =  6/a  of  (2.3). 
Specifically,  produce  the  characteristic  through  the  midpoint  of  the  west  face  of  the  control  volume, 
and  find  the  value  of  u  on  it  by  linear  interpolation  (or  extrapolation)  in  the  interval  34.  Multiply 
this  by  a  to  get  /_.  Similarly  the  characteristic  through  the  midpoint  of  the  south  face  carries  a  value 
found  from  the  data  45,  that  gives  g-  when  multiplied  by  6.  The  additional  terms  found  in  these 
fluxes  when  compared  with  the  U  scheme  can  be  regarded  as  antidiffusive  terms  which  compensate 
for  the  diffusivity  of  the  upwind  scheme.  It  is  natural  to  think  of  creating  a  high-resolution  scheme 
by  applying  limiters  to  these  terms. 

The  update  formula  for  this  scheme  is 

=  K  -  -  Us)  -  ^{U3  -  U4)  +  ^{us  -  Ui)  -  ^{uq  -  U3)] 

(3.7)  =  Uq  -  ^[^-^{U3  -  U4  -  Uo  +  Us) +  b{uo- Us)  +  a{us  -  U4)] 

(and  its  stencil  is  presented  in  Fig. 3(a). 

This  scheme  is  not  of  the  positive  type  -  the  coefficient  multiplying  U3  in  the  update  formula 
is  negative  in  our  representative  case  (3.3).  Another  defect  of  this  scheme  is  that  it  will  switch 
discontinuously  when  o  or  6  changes  sign. 
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A  steady-state  solution  for  the  homogeneous  advection  equation  obtained  by  means  of  a  zero 
cross  diffusion  scheme  will  be  second  order  accurate  [14], [5]. 

It  must  be  noted  that  there  is  an  arbitrariness  in  the  flux  formulae.  A  term  k{u3  —  U4)  can 
be  added  to  /_,  and  a  term  k{u4  —  U5)  to  5_,  without  affecting  the  update  formula,  if  k  is  any 
constant.  In  fact,  one  of  the  flux  values  is  arbitrary.  In  the  remainder  of  this  paper,  we  adopt,  as 
in  [13],  a  convention  that  if  h  >  a  then  g  is  always  given  by  the  characteristic  interpolation  (3.6), 
but  if  6  <  a  it  is  /  that  is  found  in  this  way. 


3.3.  N  scheme.  A  positive,  though  first  order  accurate,  modification  of  the  scheme  given  by  (3.6) 
is  the  following 

,  .  f!^  =  au3- lmm{a,b){u3- U4) 

^  =  bu5  —  ^min{a,b){u5  —  U4) 

or,  because  of  our  convention  (3.3), 

o'!  f-  =  -  |o(w3  -  U4) 

’  9-  =  ^^5  -  “  ’“4) 

The  corresponding  update  formula  is  (Fig.2(b)  presents  the  corresponding  stencil) 

_  .  At.,,  ,  .  , 


(3.10) 


,,n+l  _  n 

“0  ~  “0 


-[6(Tio  -  U5)  +  a(^5  -  U4)] 


This  scheme  was  introduced  by  Rice  and  Schnipke  [10]  and  was  called  the  N  (narrow)  scheme  in 
[14], [15]  for  its  narrow  stencil.  Analysis  in  [13]  shows  that  the  N  scheme  is  the  optimal  (in  terms 
of  cross  diffusion)  among  the  linear  positive  schemes  in  two  dimensions. 

It  is  interesting  to  note  that  Raithby’s  scheme  [9]  is  in  fact  identical  to  the  2D  scheme  for  the 
case  6/2  <  a  <  b.  However,  for  0  <  a  <  6/2  Raithby’s  scheme  would  correspond  to  a  blend  of  the 
2D  scheme  with  the  N  scheme.  The  weight  of  the  N  scheme  gradually  increases  with  decreasing  of 
a/6  and  becomes  equal  to  1  when  a  =  0. 


3.4.  S  scheme.  In  order  to  combine  the  positivity  property  with  second  order  accuracy  a  non¬ 
linear  scheme  has  to  be  constructed. 

A  positive  zero  cross  diffusion  scheme  can  be  constructed  in  the  way  similar  to  TVD  schemes 
in  one  dimensions.  A  “limited”  zero  cross  diffusion  correction  is  added  to  a  first  order  accurate 
positive  scheme.  Such  a  scheme  was  developed  in  [14], [15]  and  was  called  the  S  scheme.  It  is  a 
nonlinear  modification  of  the  scheme  (3.6)  and  is  given  by  the  following  fluxes 

-^^{RB43){b-a){u3-U4) 

.s  _ 
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of  a  homogeneous  advection  equation  (i.e.  it  has  a  zero  cross  diffusion)  if  is  Lipschitz  continuous 
and 


(3.18) 


^'(1)  =  1. 


Remark  3.1.  The  S  scheme  relies  on  a  5-point  stencil  (see  Fig.4(a))  -  just  one  point  more  than  the 
linear  zero  cross-diffusion  scheme  (3.6).  This  one  extra  point  was  needed  to  achieve  positivity. 

It  is  remarkable  that  the  choice  of  limiter  function  in  order  to  obtain  a  positive  scheme  with 
zero  cross-diffusion  is  dictated  by  relations  (3.17)  and  (3.18)  which  are  identical  to  those  arising 
when  constructing  one-dimensional  TVD-type  schemes.  Therefore,  most  of  the  limiters  used  in  one 
dimension  (see  [18])  can  be  used  in  the  present  context  as  well. 

Remark  3.2.  Another  zero  cross  diffusion  scheme  was  presented  by  Koren  [7] 


(3.19) 


T  N 

gt  =  9- 


The  update  formula  is 


(3.20) 


+b{uo  -  Us) 


The  advantage  of  this  scheme  is  that  it  switches  continuously  when  a  changes  sign. 

It  can  be  easily  seen  that  this  scheme  is  identical  to  the  one-dimensional  Lax-Wendroff  scheme, 
assuming  that  y  is  the  time-like  direction. 

A  detailed  analysis  of  the  class  of  zero  cross  diffusion  schemes  was  done  by  Hirsch  [5].  It  is  also 
proved  there  that  a  linear  scheme  with  zero  cross  diffusion  cannot  be  of  the  positive  type,  thus 
generalizing  Godunov’s  theorem  to  two  dimensions. 

Remark  3.3.  An  alternative  route  to  a  positive  zero  cross  diffusion  scheme  is  through  the  nonlinear 
modification  of  the  scheme  (3.19) 


(3.21) 


f-=  g  +  iX, 

£=  3^- 


However,  it  is  easy  to  see  using  (3.15)  and  assuming  a  symmetry  property  for  the  limiter,  i.e. 


(3.22) 


(and  most  of  the  commonly  used  limiters  are  symmetric,  see  [18]),  that  the  scheme  (3.21)  is  identical 
to  the  S  scheme. 

Remark  3.4.  The  choice  of  the  ratio  that  the  limiter  function  can  rely  on  is  non-unique.  For 
instance,  the  upwind  positive  scheme  with  zero  cross-diffusion  presented  by  Hirsch  &  Van  Ransbeeck 
[6]  has  a  limiter  function  with  the  following  ratio  as  its  argument 


(3.23) 


Ro34  —  - 


a(uo  -  ^3) 

b(u3  -  U4) 
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3.5.  S*  scheme.  One  of  the  main  objectives  of  this  work  is  to  establish  links  between  the  control 
volume  and  fluctuation-splitting  approaches  towards  the  construction  of  the  truly  two-dimensional 
advection  schemes.  It  is  important  for  this  purpose  to  introduce  the  scheme  S*  which  has  a  more 
‘triangular  flavor’.  The  limiter  schemes  above  are  based  on  ratios  Rijk  of  gradients  in  orthogonal 
directions.  Now  we  introduce  the  ratio  defined  as 


(3.24) 


—a{ui  —  Uj) 

(6-  a){uk  -  Uj) 


where  the  pair  of  points  j,k  are  adjacent  along  a  diagonal.  The  fraction  has  been  normalized  so 
that  we  still  have  R*jj^  ~  1  for  smooth,  nearly  steady  solutions. 

The  S*  scheme  can  be  defined  by  the  following  numerical  fluxes 


(3.25) 


/5-|^(i?S43)(6-«)(^3-«4) 


(3.26) 


a{uo  —  U4) 


{b-a)iu^-U4)'' 

The  update  formula  in  this  case  is  (the  stencil  is  presented  in  Fig.4(b)) 

r  b{uo  -  U5)  +  0(1X5  -  U4) 


(3.27) 


,,n+l  _  n 
“0  ~  ^0 


4-^^(jR043)(^-  o)(«3  -  U4) 
“2^(-R75o)(^  -  a)(ixo  -  U5)] 


Note  that  the  following  identity  holds 


(3.28) 


'i{Ro43){b  -  a)(iX3  -  1x4)  =  -^^^0(1x0  -  U4). 

no43 


Lemma  3.5.  If  the  limiter  4'  =  '^{R)  satisfies  the  following  inequality 

(o  on^  n  /  ^  O  O  /  iTf  X  ^  O 


(3.29)  0  < 

then  the  S*  scheme  is  of  positive  type. 


<  2,0  <  $(i?)  <  2. 


Proof.  Using  the  identity  (3.28)  the  update  formula  for  the  S*  scheme  can  be  rewritten  as  follows 

.  r  {b-  a)  (1  -  |$(i?75o))  ("fio  -  U5)  ■ 

(3.30)  ■ 


..n+l  _  ..n  _ 

Uq  —  ^0 


It  is  obvious  from  (3.30)  that  the  scheme  is  positive  if  the  inequality  (3.29)  holds.  □ 

Lemma  3.6.  7/4' (i?)  is  Lipschitz  continuous  and 
(3.31)  ^(1)  =  1, 

then  the  S*  scheme  has  a  zero  cross  diffusion. 

Proof.  This  lemma  is  a  simple  corollary  of  Lemma  A.l  concerning  the  more  general  scheme  con¬ 
structed  for  the  inhomogeneous  equation.  □ 
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Figure  5.  Illustration  for  construction  of  fluctuation-splitting  schemes. 


4.  Fluctuation  splitting  approach 


In  this  section  we  present  a  brief  description  of  a  fluctuation  splitting  approach.  For  more  details 
see  [3], [12], [17].  The  schemes  are  designed  for  use  on  unstructured  grids  and  their  description  is 
initially  given  in  that  context.  Later,  to  compare  with  the  control-volume  approach,  we  specialize 
to  structured  grids,  and  later  still  we  reconsider  the  general  case. 

Consider  a  numerical  solution  of  the  two  dimensional  linear  constant  coefiicient  advection  equa- 


(4.1)  ut  +  A-Vu  =  s. 

The  grid  is  taken  to  be  an  arbitrary  triangulation  of  the  domain.  A  typical  element  T  of  such  a 
grid  is  shown  on  Fig. 5.  The  integral  of  Uf  over  the  triangle  T  is  called  the  fluctuation 


Utdxdy. 


Taking  into  account  Eq.(4.1)  and  applying  Gauss’  theorem  we  obtain 

(4.3)  (f)^  =  ^  uX-  dn-  JJ  s  dxdy, 

where  dT  is  the  boundary  of  the  triangle  T  and  dn  is  the  inward  scaled  normal  to  an  element  of 
the  boundary.  Assuming  that  the  source  term  s  varies  linearly  within  triangle  T  we  get 

(4.4)  <7T  =  JJ  s  dxdy  = 

where  St  is  the  area  of  the  triangle  T.  We  shall  return  to  the  inhomogeneous  case  in  Appendix  B. 
Here  we  assume  that  s  =  0  (and  therefore  ax  =  0). 

(4.5)  Ut  X  ■  Vu  =  0. 

Assuming  also  that  u  varies  linearly  within  triangle  T  we  get 

1  ^  1  --.1 

(4.6)  4>^  =  2(^1  +  ^2)A  •  ^3  +  2(^2  +  «3)A  •  ni  +  -(ui  +  ^3) A  •  n2, 
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Figure  6.  Median  dual  cell. 


where 

(4.7) 

and  rii  is  the  inward  normal  to  the  side  Ei  scaled  with  its  length.  Note  that 


(4.8) 

and  therefore 

(4.9) 

This  allows  (4.6)  to  be  rearranged  as 

(4.10) 


=  0 


i=l 


— 0- 


i=l 


4>'^  =  -'^kiUi. 


i-\ 


Several  alternative  expressions  can  be  obtained  for  (j)^  using  (4.9),  for  example, 

<f)^  =  -k2{u2  -  Ui)  -  k3{u3  -  Ui) 

=  -k3{u3  -  U2)  -  ki{ui  -  U2) 

(4.11)  =  -ki{ui  -  U3)  -  k2{u2  -  U3) 

The  computed  fluctuation  should  then  be  distributed  (split)  among  the  vertices  of  the  triangle  T 

(4.12)  SiUi  :=  SiUi  +  i  =  1,  2,  3, 

where  5,-  is  one-third  the  area  of  all  the  triangles  having  as  a  vertex  (the  area  of  so-called  median 
dual  cell  of  area  Si  around  i,  see  Fig. 6)  and 

(4.13)  = 

i=l 

The  latter  is  necessary  for  the  discretization  to  be  conservative  (see  [3]). 

After  adding  contributions  from  all  the  triangles  T  having  a  common  vertex  at  point  i,  we  obtain 
the  following  scheme 


(4.14) 
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(a)  Type  I  triangle  (b)  Type  II  triangle 

Figure  7.  Two  possible  triangle  orientation  with  respect  to  the  flow  direction. 


We  assume  here  that  each  triangle  sends  contributions  to  its  own  vertices  only.  Therefore,  the 
resulting  scheme  for  Ui  has  a  compact  stencil  restricted  to  at  most  the  immediate  neighbours  of  ui. 
The  remaining  question  is  how  to  chose  aj . 

An  additional  consideration  is  positivity  of  the  constructed  scheme.  We  will  define  a  more 
restrictive  condition  of  local  positivity  that  requires  contributions  of  each  triangle  separately  to 
be  positive. 

Definition  4.1.  In  equation(4.12)  choose  from  the  preceding  list  the  definition  of  <j)^ ,  so  that 

(4.15)  Si'tii  I—  S-iiii  (X-I^kj(uj  r/j)  -h 

Then  the  scheme  is  said  to  be  locally  positive  if  the  three  quantities 

-OiikjAt,  —OikkAt,  1  +  ai{kj  +  kk)At 

are  all  positive 

This  condition  is  easier  to  implement  because  it  retains  the  basic  property  of  a  fluctuation¬ 
splitting  scheme  that  all  triangles  can  be  processed  without  reference  to  any  other  data.  However, 
it  is  more  restrictive  than  necessary  because  it  does  not  recognise  that  compensating  changes  from 
other  triangles  may  restore  positivity.  It  is  the  condition  currently  used  in  design  of  fluctuation 
splitting  schemes.  Application  of  the  non-local  positivity  condition  may  allow  addition  of  artificial 
compression  to  the  scheme.  This  will  be  discussed  in  §5.2. 

Definition  4.2.  The  fluctuation-splitting  scheme  is  called  linearity  preserving  if  whenever  the  fluc¬ 
tuation  on  the  triangle  T  vanishes  then  the  scheme  leads  to  a  zero  update  in  each  of  the  three  vertices 
of  the  triangle. 

It  was  observed  by  Deconinck  et  al  [3]  that  linearity-preserving  schemes  give  solutions  that  are 
second-order  accurate  in  the  steady  state. 

4.1.  N  scheme.  This  linear,  positive,  fluctuation  splitting  scheme  is  called  the  N  scheme  because 
it  was  found  by  Roe  [12]  that  it  is  closely  related  to  the  control  volume  N  scheme.  (We  shall  return 
to  this  point  in  §5.1.1).  Note  that  there  exist  two  possibilities  regarding  ki  for  i=  1,2,3  (because 
of  Eq.(4.9)):  either  one  of  them  is  positive  and  the  rest  are  negative  or  two  of  them  are  positive 
and  one  negative.  These  two  situations  correspond  to  a  triangle  with  either  one  inflow  face  or  two 
inflow  faces.  Without  loss  of  generality  we  can  consider  the  following  two  cases  (see  Fig. 7) 

(a)  :  ki  <  Q,  k2  <  0,  fcs  >  0, 

(b)  :  ki  >  0,  k2  >  0,  ks  <  0. 
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The  obvious  choice  for  the  first  case  (Type  I  triangle)  is  to  send  the  entire  fluctuation  to  the 
downstream  vertex  3  (i.e.  setting  =  aj  =  0,  0:3  =  1  ) 

(4.16)  =  S3U3  -  At[ki(ui  -  U3)  +  ^2(^2  -  U3)] 

Although  interesting  schemes  can  be  created  that  do  not  follow  this  rule,  we  will  assume  for  the 
remainder  of  this  paper  that  all  Type  I  triangles  are  dealt  with  like  this. 

Case  (b)  is  more  complicated.  If  the  fluctuation  is  distributed  according  to  the  following  formulae 

-  Af[fci(«i  -  U3)] 

^  ^  S2U2'^^  =  52«2  -  At[A:2(ti2  -  Ms)] 

the  resulting  scheme  is  obviously  positive  and  this  will  define  the  N  sceme.  Strong  convergence  of 
the  N  scheme  on  general  triangulations  is  proved  by  Perthame  [8].  However,  it  does  not  have  the 
LP  property.  This  is  because  the  contributions  to  the  residuals  at  both  vertices  1  and  2  may  be 
different  from  zero  even  if  the  fluctuation  on  triangle  T  vanishes. 

4.2.  NN  scheme.  It  was  pointed  out  in  [11], [3], [12]  that  replacing  the  advection  velocity  A  by 
its  normal  to  the  level  lines  component  A„  normal  to  the  level  lines  (or  tangential  to  the  gradient 
direction)  of  u  within  triangle  T 


(4.18) 


(A  •  Vu) 


does  not  affect  the  residual.  It  is  also  obvious  that  substituting  the  advection  velocity  by  the 
following 

(4.19)  A*  —  A„  +  9X 

does  not  affect  the  residual  either.  Defining 


(4.20) 


k-  =  2^* 


5im^+i  =  Siu^  -  At[k*{ui  -  U3)] 

^  ‘  '  52M2'*'^  =  52 -  At[k2{u2  -  M3)] 

one  can  observe  that  if  0  =  0,  the  NN  scheme  (4.21)  will  be  linearity  preserving.  This  is  because 
A„  vanishes  at  the  steady  state  (and  therefore  fcj  and  vanish  in  this  case).  However,  vector  A„ 
may  extend  out  of  triangle  T  (Fig. 8).  Therefore,  the  scheme  in  this  case  is  not  positive,  i.e.  one  of 


the  it’s  in  (4.21)  will  be  negative.  To  obtain  positivity  of  the  scheme  in  this  case  d  should  be  taken 
such  that  the  vector  A*  will  be  directed  along  the  face  of  triangle  T.  This  means  that  the  entire 
residual  should  be  sent  to  only  one  of  the  outflow  nodes. 

Because  the  effective  advection  speed  is  data  dependent  (when  0  =  0  it  is  the  drift  velocity  of 
the  contours)  the  NN  scheme  is  nonlinear  in  a  way  that  permits  it  to  escape  Godunov’s  Theorem. 
Although  the  mechanism  involved  appears  very  different  from  that  employed  by  limiter  schemes, 
we  will  see  in  the  next  section  that  a  unification  is  possible. 

It  was  pointed  out  in  [3]  that  when  the  solution  is  almost  a  constant,  the  expression  for  A„  (4.18) 
becomes  ill-defined,  which  may  affect  the  convergence.  A  way  to  overcome  this  difficulty  suggested 
in  [3]  was  to  switch  gradually  from  A*  to  A  in  regions  where  the  gradients  are  small. 


4.3.  Level  scheme.  It  was  argued  by  Roe  in  [12]  that  the  NN  scheme  is  not  optimal  in  terms  of 
the  maximal  allowed  time  step  and  another  scheme  for  Type  II  triangles  was  suggested  there.  This 
scheme  was  later  named  the  Level  scheme. 

A  brief  description  of  the  algorithm  is  given  here.  The  reader  is  referred  to  [12]  for  the  complete 
derivation. 

The  residual  is  distributed  among  the  nodes  according  to  the  following  formulae: 

-!-  Atai(f>T 

^  =  S2U^  +  Ato2<^T 

It  is  very  interesting  for  the  discussion  in  §5.1.3  to  note  that  in  [12]  the  following  quantity 


(4.23) 


Ui  -  U3 
«2  -  '^3 


was  suggested  to  distinguish  between  two  cases 

(a)  :  A  >  0,  i.e.  the  level  line  passes  through  triangle  T, 

(b)  :  A  <  0,  i.e.  the  level  line  lies  outside  of  triangle  T. 

In  case  (a)  only  one  vertex  is  updated:  the  vertex  that  the  characteristic  line  passes  closer  to  it 
than  the  level  line. 

In  case  (b)  both  vertices  are  updated.  The  residual  is  distributed  according  to  the  following 
formulae 


(4.24) 


1  5l(til-«3)+52(«2-«3) 

_ 

2  5l(^il-1i3)^-S2(^i2-^^3) 


This  choice  of  a’s  is  argued  in  [12]  to  be  optimal  in  terms  of  the  maximal  allowed  time  step.  Also, 
this  choice  preserves  the  direction  of  the  level  lines;  hence  the  name  of  the  scheme. 

However,  Deconinck  [1]  has  observed  that  the  Level  scheme  may  switch  discontinuously  if  A  is 
almost  parallel  to  one  side  of  the  triangle,  say  side  1,  but  (^2  — 'U3)/('ni  —us)  >  0.  This  demonstrates 
that  A  was  an  unhappy  choice  of  parameter. 


5.  Unified  approach 

It  is  shown  in  this  section  that  there  is  more  in  common  between  the  two  previously  reviewed 
approaches  than  just  a  philosophy  to  make  constructive  use  of  nonlinearity.  In  fact  the  fluctuation 
splitting  methodology  and  control  volume  compact  scheme  methodologies  can  be  seen  as  two  dual 
ways  to  interpret  the  same  approach. 

Consider  the  following  scheme  on  Type  II  triangles 

.  ^  =  Siu^  +  At[k-i{ui  -  U3) +  '^{Ri32)k2{u2  -  U3)] 

^  '  S'2M2'^^  =  S2U'^  +  At[k2{u2  -  U3)  -  fl'(i?132)fc2(w2  “  “s)] 
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where  Rijk  now  compares  contributions  to  the  fluctuation  from  two  sides  of  the  triangle,  i.e. 

„  kdui  —  u-i) 

(5.2)  Rijk  =  -ttH - 4- 

Again,  we  have  Rijk  —  1  for  smooth,  nearly  steady  flow. 

^  kiiui  —  Us) 

=  ~hU-us) 

The  conservation  property  of  this  scheme  is  obvious  from  (5.1). 

We  have  the  identity 

Cr  iTfCD  \  7.  ^  —  ^(.^132)  \ 


^132  —  —  ■ 


^(.^132) ^2 ('^2  -  Us)  = 


-k-i{ui  -  Us) 


Theorem  5.1.  If  the  limiter  'i’{R)  satisfies  the  following  inequality 

(5.5)  0  <  <1,  0  <  ^{R)  <  1, 

then  the  NNL  scheme  is  (locally)  positive. 

Proof  Using  the  identity  (5.4)  the  scheme  (5.1)  can  be  rewritten  in  the  following  form 

Siu^+^  =  5i<  +  At[l  -  -  us) 

^  ’  S2U^+^  =  S2U^  +  Af[l  -  '^{RlS2)]k2{U2  -  Us) 

It  can  be  concluded  from  the  latter  that  the  NNL  scheme  is  positive  if  inequality  (5.5)  holds.  □ 

Remark  5.2.  By  contrast  to  the  control  volume  schemes,  it  is  here  possible  only  to  use  limiters  for 
which  the  upper  bound  on  is  equal  to  1.0  rather  than  2.0.  That  is,  the  limiter  cannot 

be  compressive.  We  show  later  that  compressive  limiters  can  be  allowed  within  the  fluctuation 
splitting  scheme  given  sufficient  information  about  mesh  connectivity. 

Theorem  5.3.  The  NNL  scheme  is  linearity  preserving  if 


(5.7) 

for  some  positive  constant  C . 
Proof.  We  have 

^Sui 

^  -  S2U^  =  5u2 


\^{R)-1\<C\R-1\. 


At[ki{ui  -  Us)  +  4'(J?i32)A:2(u2  -  Us)]  =  p  +  '^{R)Q 
At[k2{u2  -  us)  +  ^^h{u^  -  us)]  =Q-  ^{R)Q 


We  wish  to  show  that  if  F  +  Q  =  0  then  Sui  =  5u2  =  0,  or  equivalently 

|5ui  —  ^■U2|  <  K\Sui  +  6u2\, 
for  some  positive  constant  K.  Now,  if 

I^'-IKCIF-Il 

then 

2|^-1||Q|  <  2C\P  +  Q\, 

|F  +  Q|  +  2|^-1||Q|  <  (2C+1)|F  +  Q|, 
|F  +  Q  +  2(^-l)Q|  <  (2C  +  1)|F  +  Q|, 
|F-g  +  2^Q|  <  (2C+1)|F  +  Q|. 


This  last  statement  translates,  with  K  =  2C  +  1,  into  the  one  we  want  to  prove. 


2  1  8 


Figure  10.  Possible  triangulation  of  the  Cartesian  grid. 


Remark  5-4.  It  is  interesting  to  note  that  the  NNL  scheme  with  the  “minmod”  limiter  is  identical 
to  the  Level  scheme  in  the  case  when  the  level  line  lies  inside  the  triangle.  To  illustrate  this  assume 
that  Ri32  >  1.  For  the  minmod  limiter  this  means  that  V['(i?i32)  =  1-  It  follows  from  (5.1)  that  the 
entire  fluctuation  contributes  to  the  solution  update  at  the  node  1  -  exactly  what  the  Level  scheme 
corresponds  to. 

5.1.  Links  to  the  control  volume  approach. 

5.1.1.  N  schemes.  Consider  the  grid  as  illustrated  in  Fig.9.  The  fluctuations  for  the  triangles 
A,  B,  C  are  given  by  the  following  expressions 

=  -^[{b  -  a){u3  -  U4)  +  a{uo  -  u^)] 

(f>^  =  -^[<^(^0  -  '“4)  +  (^  -  q)(^o  -  ^s)] 

(t>^  =  -  -  U5)  +  {b-  a){uo  -  Us)] 

The  underlined  terms  here  contribute  to  the  residual  at  the  node  0  according  to  (4. 16), (4. 17). 
Assembling  all  the  contributions  the  following  update  is  obtained 

(5.10)  +  +  a(u4  -  Wo)] 

It  was  pointed  out  by  Roe  [12]  that  (5.10)  is  identical  in  this  case  to  the  update  formula  for  the 
control  volume  type  N  scheme  (3.10).  That  is  why  these  two  schemes  are  called  N  schemes.  Another 
interesting  observation  made  in  [12]  is  that  if  the  N  scheme  is  constructed  on  the  triangulation 
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illustrated  on  Fig.  10 


=  -|[(6  -  a)(w3  -  Ui)  +  a{uo  -  -U4)] 

=  -•|[q(^o  -  '^4)  +  (fe  -  Q)(txo  -  Us)] 

=  -2[“(^7  -  %)  +  ib-  a){uo  -  -Us)] 
it  can  be  seen  from  the  update  formula 

(5.11)  ^0+^  "US'  +  -  "Mo)  +  a{u3  -  -uo)] 

that  the  scheme  in  this  case  is  identical  to  the  dimensional  upwinding  (3.5). 

5.1.2.  2D  scheme.  Consider  again  triangulation  illustrated  on  Fig. 9.  The  fluctuations  on  the  tri¬ 
angles  A,  B,  C  are  again  given  as  follows 

=  -f  [(&  -  ffl)(M3  -  ^4)  +  ^(^0  -  •^4)] 

=  -t[Q(tro  -  U4)  +  {b-  a){uo  -  -U5)] 

(ff  =  -|[a(ti7  -  U5)  +  {b-  a){uo  -  ws)] 

Note,  that  now  only  triangles  A,  B  participate  in  the  molecule  at  node  0,  however  they  contribute 
their  entire  fluctuations.  This  would  correspond,  in  the  fluctuation-splitting  approach,  to  a  rule 
that  each  triangle  only  updates  it’s  most  downwind  node  (the  one  for  which  A  •  x  is  maximum). 
The  update  formula  in  this  case 

.  ,  „  AtJa-b).  .  (b—a),  .  (a  +  b). 

(5.12)  rio'*'  =  <  +  -j-i  2  ~  '  2  ^  ^  2  ~ 

is  identical  to  the  2D  control  volume  scheme  (3.7). 


5.1.3.  Linearity  preserving  positive  schemes.  Consider  triangles  A,B,C  as  illustrated  in  Fig. 9.  De¬ 
fine 


(5.13) 


_  a{uo  -  Uj) 

■^043  —  /I  \  /  \  )  ^750  — 

(6  -  a)  [Us  -  U4) 


a{u7  -  -Us) 

[b  -  a){uo  -  U5) ' 


The  fluctuations  on  these  triangles  are 


<f)^  =  -f  {[1  -  ^^(i2o43)](^»  -  a){'U'3  -  U4)  +  a{uo  -  U4)  +  fl>(j?043)(6  -  a)(tr3  -  ^4)} 

(5.14)  {a(uo  -  ^4)  +  {b-  a){uo  -  -us)} 

=  “2  {“('“7  -  '“5)  +  4’(i?75o)(^'  -  a)(Mo  -  Ms)  +  [l  -  ( ^750)] (6  -  a)  (rtp  -  Ms)} 


The  underlined  terms  contribute  to  the  residual  at  the  point  0 


r  {{b  -a)-  |^'(i?75o)(^»  -  a))  (m5  -  Mo) 
(5.15)  Uq'^'^=Uq  +  —  -f  (a-b  |^'(fio43)(^»  -  a))  (m4  -  Mo) 

[  +  (-|^(i?043)(&  -  a))  {U3  -  Uo)  _ 


which  is  identical  to  the  update  formula  for  the  S*  scheme  (3.27)! 


5.2.  Artificial  compression.  Note  that  it  was  possible  to  use  a  compressive  limiter  ^  for  the 
S*  scheme.  This  was  possible  to  show  by  requiring  non-local  positivity,  which  is  a  less  restrictive 
requirement.  If  we  were  to  restrict  our  attention  only  to  structured  grids,  we  could  use  the  equiva¬ 
lence  of  the  S*  and  NNL  schemes  to  justify  using  compressive  limiters  in  the  NNL  scheme  also. 
We  now  discuss  how  compressive  limiters  may  be  employed  when  we  know  something  about  the 
connectivity  of  the  grid. 

Consider  triangle  C  (Fig. 9) .  Collect  from  triangles  B  and  C  those  terms  in  the  fluctuations  that 
contribute  to  the  update  at  node  0  and  contain  a  factor  {u^  —  "Uo)- 


compressive  limiter 


.  1.  -  J  _ ] 
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(5.16) 


Figure  11.  General  triangulation. 


4>q+4>0  =  ~  ^(^75o)](^- a)(rf5  -  +  (^>  -  a)('M5  -  Wo)} 


Similarly,  collect  from  triangles  D  and  C  those  terms  in  the  fluctuations  that  contribute  to  the 
update  at  node  7  and  contain  a  factor  (ws  —  uj). 


4>?  +  (t>?  =  |{[1  -  -  a){u5  -  uj)  +  {b-  a){u5  -  uj)} 

^  /t7R0 


(5.17) 


R750 

'^iR75o) 


](6-  a)(u5  -  U7) 


It  is  obvious  that  both  these  contributions  will  be  positive  if  the  following  inequality  holds 

/cr  1  o>  f\  ^  /  o  n  ^  iTr/"  r?^  /  o 


(5.18) 


<2,  0  <  ^{R)  <  2. 


This  means  that  a  compressive  limiter  can  be  used  in  this  case.  However,  this  case  is  particularly 
simple:  the  grid  is  nicely  structured  and  the  advection  velocity  is  constant.  The  case  of  unstructured 
grids  and/or  variable  advection  velocity  is  more  complicated.  (  The  nonlinear  case  can  be  reduced 
for  this  purpose  to  the  variable  velocity  case  by  using  the  multidimensional  conservative  linearization 
(see  [2])). 

Consider  this  general  situation  as  illustrated  in  Fig.ll.  Assume  that  triangle  C  is  of  Type  II, 
while  the  neighboring  (through  the  inflow  faces  of  C)  triangles  B  and  C  are  of  Type  I.  Collect,  as 
before,  the  contributions  to  the  updates  at  nodes  0  and  7  containing  factors  (its  —  wq)  and  (its  —  W7). 


(5.19) 


(t>o+<t’o  =  -  Wo)  +  A:^(its  -  Ito)} 

=  ^[{^0  +  ^0  )  -  ^(^75o)A:^](w5  -  Wo) 


(5.20) 


4'7+(f>?  =  |{[1-^-^^]^?(w5-W7)  +  ^?(ws-  W7)} 

^  ^750 

^  tt7so 


<  Tac,  0  <  W(i?)  <  7ac. 


0  14^  1  X,  R 

Figure  12.  Positivity  region  for  the  limiter  used  in  the  NNL  scheme. 

where  fc’s  are  defined  according  to  (4.7)  (superscripts  indicate  the  triangle).  Defining  the  following 
‘artificial  compression  factor’ 

kC  xuD  I.C  \  uB 

(5.21)  7ac  =  niin(  ) 

we  can  state 

Theorem  5.5.  The  NNL  scheme  is  of  the  positive  type  if  the  limiter  '^{R)  employed  in  triangle  C 
satisfies  the  following  inequality 

(5.22)  0  <  <  7ac,  0  <  W(i?)  <  7ac. 

Proof.  The  proof  is  obvious  from  the  previous  argument.  □ 

The  positivity  region  (5.22)  for  the  limiter  function  used  with  the  NNL  scheme  is  illustrated  on 
Fig. 12.  A  class  of  limiters  defined  by 

(5.23)  =  max(0,  mm{<j)R,  1),  min(i?,  4>)) 
with  the  parameter  ^  satisfying 

(5.24)  l<<t><  lac- 
appears  to  suit  this  case  very  well.  It  is  obvious  that  the  choice 

(5.25)  =  1 

gives  the  minmod  limiter  which  corresponds  to  the  lower  bound  of  the  positivity  region  (Fig.l2). 
The  choice 

(5.26)  (j^  —  lact 

results  in  the  “generalized”  highly  compressive  Roe  superbee  limiter  and  corresponds  to  the  upper 
bound  of  the  positivity  region  (see  Fig. 12). 

The  amount  of  artificial  compression  that  the  scheme  is  allowed  to  use  should  probably  be 
problem  dependent.  Usually  it  does  not  make  much  difference  for  the  resolution  of  strong  shocks 
and  may  lead  to  some  undesirable  effects  in  the  smooth  regions  (for  example,  an  advected  sinusoidal 
wave  may  start  looking  like  a  square  wave).  However,  it  can  make  a  great  difference  in  resolving 
contact  discontinuities. 
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Figure  13.  Regular  upwinding  and  N  scheme. 

6.  Some  numerical  experiments 

In  this  section  we  present  some  numerical  experiments  which  illustrate  the  performance  of  the 
genuinely  two-dimensional  advection  schemes.  Two  important  aspects  are  considered:  resolution  of 
discontinuities  and  the  accuracy  in  smooth  regions.  All  testcases  presented  in  this  section  consider 
the  linear  advection  equation 

=  0. 

6.1.  Resolution  of  discontinuities.  An  extensive  set  of  numerical  experiments  with  the  5  scheme 
concerning  the  resolution  of  discontinuities  is  presented  in  [15].  Here  we  shall  illustrate  this  ability 
of  the  S*  (or  NNL)  in  comparison  to  the  S  scheme. 

We  consider  here  the  advection  equation  (6)  on  the  square  [0, 1]  X  [0, 1].  The  boundary  conditions 
and  the  exact  solution  steady-state  solution  are  given  by 

u  =  H{y-  (2x-2)), 

where  H  is  the  Heaviside  function.  The  grid  used  is  50  X  50  points. 

Figuresl3-16  present  contour  plots  of  the  numerical  solution  to  this  problem  obtained  using 
different  schemes.  The  performance  of  the  first  order  regular  upwind  and  the  N  schemes  is  illustrated 
by  Fig.l3.  It  can  be  observed  that  the  N  scheme  is  superior  to  the  dimensional  upwinding  in  terms 
of  discontinuity  resolution.  However,  the  width  of  the  numerical  layer  representing  the  discontinuity 
is  in  both  cases. 

Fig. 14  presents  the  results  obtained  by  the  S  and  the  S*  (which  is  the  same  as  NNL  scheme  in 
this  case)  employing  the  non-compressive  mmmod  limiter.  Fig. 15  corresponds  to  the  same  schemes 
except  that  the  van  Leer  (harmonic  mean)  limiter  is  used.  Experiments  with  the  highly  compressive 
superbee  limiter  are  presented  in  Fig. 16.  It  can  be  seen  that  the  discontinuities  are  represented  by 
sharp  layers  with  no  oscillations  in  the  numerical  obtained  by  both  the  S  and  the  S*  schemes  with 
all  the  three  different  limiters.  Another  conclusion  is  that  the  artificial  compression  may  have  a 
strong  effect  on  the  sharpness  of  such  a  layer.  It  can  also  be  observed  that  the  S*  is  slightly  superior 
to  the  S  scheme  in  resolving  contact  discontinuities.  This  can  be  explained  by  its  slightly  narrower 
stencil. 
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Figure  14.  S  and  NNL  scheme  with  “minmod”  limiter. 


Figure  15.  S  and  NNL  scheme  with  van  Leer  limiter. 

6.2.  Accuracy  in  smooth  regions.  The  purpose  of  the  experiments  presented  here  is  to  verify  the 
accuracy  of  the  S*  (or  the  NNL  scheme  on  structured  grid)  in  the  smooth  regions  and  to  compare 
it  to  that  of  the  S  scheme.  An  extensive  set  of  numerical  experiments  with  the  S  employing 
different  limiters  was  presented  in  [15].  Here  we  restrict  our  attention  to  the  minmod  limiter  since 
the  S*  scheme  (as  well  as  the  NNL  scheme  in  general  unstructured  grid  case)  satisfies  the  local 
positivity  property  in  this  case.  Therefore,  the  treatment  of  each  triangle  is  purely  local.  This 
results  in  a  significantly  more  economical  numerical  algorithm,  though  at  the  expense  of  slightly 
inferior  resolution  of  contact  discontinuities. 

The  testcase  considered  is  again  the  advection  equation  (6)  on  the  square  [0,1]  X  [0,1].  The 
boundary  conditions  (and  the  exact  steady-state  solution)  are  given  by 

u  =  sin  (4a;  —  2y) 
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Figure  16.  S  and  NNL  scheme  with  Roe  superbee  limiter. 


gridsize 

12x12 

25x25 

50x50 


S  scheme 
1.40  X  10-^ 
3.38  X  10-^ 
8.54  X  10-^ 


S*  scheme 
1.17  X  10-^ 
2.69  X  10-^ 
6.73  X  10-4 


Table  1.  L\  solution  error  norm  (measured  on  the  [0,0.75]  X  [0,0.75]  subdomain  to 
avoid  the  influence  of  numerical  boundary  layers)  for  the  problem  Ut-\-Uy  +  .bUx  =  0 
with  the  boundary  conditions  (and  the  exact  steady-state  solution)  u  =  sin(4x  — 2t/). 


The  Li  solution  error  norm  measured  on  the  subdomain  [0,  0.75]  X  [0,  0.75]  (to  avoid  any  influence 
of  the  numerical  boundary  layers)  is  presented  on  Table  1  for  the  cases  of  the  S  and  S*  schemes 
and  three  different  levels  of  resolution. 

The  experiments  with  both  schemes  demonstrate  second  order  accuracy  at  the  steady-state  of 
the  both  schemes.  The  solution  error  in  the  case  of  the  S*  scheme  is  slightly  smaller  than  for  the 
5  scheme.  The  explanation  for  this  is  the  slightly  narrower  stencil  of  the  5*  scheme. 

7.  Discussion  and  conclusions 

This  work  establishes  a  link  between  two  types  of  genuinely  two-dimensional  advection  schemes. 
The  truly  two-dimensional  schemes  of  the  control  volume  type  were  constructed  in  [14], [15].  The 
fluctuation-splitting  type  schemes  were  introduced  in  [11], [3], [12], [17].  The  S*  scheme  presented 
in  this  work  can  be  interpreted  as  a  representative  of  both  classes.  This  allowed  us  to  transfer 
ideas  and  insights  between  the  two  approaches  and  to  combine  the  desirable  properties  of  both 
classes  in  a  single  scheme.  This  is  accomplished  by  the  NNL  scheme,  which  on  one  hand  (as  a 
fluctuation-splitting  scheme)  is  formulated  on  unstructured  triangular  grids.  Thus,  it  allows  for  a 
simple  treatment  of  complex  geometries.  On  the  other  hand  the  NNL  scheme  inherits  the  following 
properties  of  the  control  volume  type  genuinely  multidimensional  advection  schemes: 

•  to  use  regular  limiter  functions  in  order  to  formulate  the  fluctuation-splitting  schemes  for 
unstructured  triangular  meshes; 

•  possibility  to  use  artificial  compression,  which  can  improve  the  resolution  of  discontinuities; 
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•  to  formulate  a  fluctuation-splitting  scheme  having  the  LP  property  for  the  advection  equation 
with  nonzero  forcing  term. 

One  of  the  important  objectives  of  this  research  direction  is  to  construct  highly  efficient  and 
robust  multigrid  steady-state  solver  for  the  compressible  flow  problems.  The  main  motivation  for 
constructing  the  multidimensional  schemes  of  the  control  volume  type  in  [14], [15]  was  to  obtain  a 
high-resolution  scheme  such  that  Gauss-Seidel  relaxation  (the  simplest  and  very  efficient  smoother) 
is  stable  when  applied  directly  to  to  the  resulting  discrete  equations  (which  is  not  the  case  for 
dimensionally-split  high-resolution  schemes  [16]).  This  made  it  possible  to  construct  a  highly  effi¬ 
cient  multigrid  steady-state  solver.  The  genuinely  two-dimensional  advection  schemes  formulated 
in  this  work  (5*  and  NNL)  inherit  this  property  from  their  control  volume  type  predecessors. 
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Appendix. 


Appendix  A.  S  scheme  -  general  inhomogeneous  case 
Consider  the  following  scheme 

,  .  /^^_  =  aus  -  l{b{u3  -  U4)  -  hsi) 

bu5  —  ^{a{u5  —  U4)  —  hs^_). 

It  has  been  shown  in  [14], [15]  that  the  2D  scheme  is  second  order  accurate  at  steady  state  for 
inhomogeneous  equation.  The  definition  of  is  non-unique  (see  [14], [15]).  We  consider  here 

one  particular  choice  which  appears  to  be  of  particular  importance  for  the  purpose  of  this  work 
(see  Appendix  B): 

.  s  f  2  1 

(A.2)  s_  —  -S3  -f  -S4 

(A. 3)  si  =  |s5  +  is4 

If  one  is  content  to  obtain  first  order  accuracy  the  previously  defined  N  scheme  can  be  used. 
However,  the  purpose  of  it  here  is  to  serve  as  an  intermediate  step  towards  the  construction  of  the 
second  order  accurate  S*  scheme.  Therefore,  we  shall  modify  it.  For  our  representative  case  a  <  b, 
putting 

(A.4)  /3  =  min(l,^)  =  ^  <  1, 


and  defining 


the  N  scheme  can  be  given  by: 


The  S*  scheme  is  defined  by 


s-  —  ^(^0  +  S3  +  S4) 


f^_  =  au3  -  Ha(u3  -  U4)  -  h[si  -  (1  -  /9)s_]} 
g^_  =  bus  —  2(®(^5  “  '“4)  — 


=  f^_  -  i«f(i?043)((^>  -  o)(«3  -  U4)  -  (1  -  P)hs  ) 


9^-  =  9^- 


where 


-R043  — 


The  update  formula  in  this  Ccise 


a(tio  —  "“4)  “  /3h.s_ 

(6  —  o)(u3  —  'M4)  —  (1  —  /3)/iS_ 


,n+l  _  ^,n 


Ug  -f-  — {  b(u5  -  Uo)  4-  a(u4  -  U5) 


-^'(i?75o)[(h  —  u)(u5  —  Uo)  +  h(l  —  /3)s- 


+-^(i?043)[(6  -  a)(u4  -  y-s)  +  h(l  -  /?)s+] 
+^[si  -  si  +  +  (1  -/3)(s+  -s_)]} 


or 


Uq'^''-  =Uq  +  —{  6(rt5  -  tio)  +  a(u4  -  fis) 

--^(/?75o)[(6  -  a){u^  -  Mo)  +  h{l  -  I3)s-] 

+^^'(i?043)[(&-  «)(m4  -  Ms)  +/r(l  -  /3)5+] 

(A. 10)  +  (^  ""  ^)'®3  +  (2  —  P)s4  +  (1  +  P)s5  +  /^^r]} 

Lemma  A.l.  If^  =  '^{R)  is  Lipschitz  continuous  and 
(A.ll)  ’J'(l)  =  l, 

then  the  S*  scheme  is  second  order  accurate. 

Proof.  Note  that 

iff  -  /f )  -  (/f  -  /!^) 

—  (1  —  ^(-R75o))((^  —  m)(mo  —  Ms)  — /3/iS+) 

-(1  -  «'(i2o43))((&-  m)(m3  -  M4)  -  Phs^) 

=  /i(^(i?043)  -  ^(-R75o))((^  -  a)  (My  -  .5hUyy)  -  13s) 

+/l^(l  -  ^'(/2o43))((^-  ofUxy  -  (Ssx) 

(A.12)  +Oih^) 

The  second  order  accuracy  of  the  S*  scheme  requires: 

(A. 13)  1  —  ^'(i?043)  =  ^(^) 

and 

(A.14)  ^(i?043)  -  ^{Rrso)  =  0{h^) 


(A.15) 

Similarly 

(A.16) 

Therefore 

(A.17) 


=  1-h 


a(uo  —  U4)  —  0hs- 

{b  -  a)(u3  -  U4)  —  (1  —  fi)hs- 
a(uo  —  M4)  +  (6  —  a)(u3  —  ■M4)  -  hs- 
(b  -  a){u3  —  U4)  —  (1  —  fi)hs- 
{aUx  +  buy  —  s)  —  .bhbuyy  —  .bhau^x  —  hbu^y  —  .bhs^ 
{b  —  a)uy  —  (1  —  S3)s  +  0{h) 

^  .bbUyy  .baUxx  bUxy  .bSx  ^  C^^h^) 


+  Oih'^) 


(b  —  a)uy  —  (1  —  I3)s 


iJrso  =  1  -  h  ■5H,  +  .5a«..  +  .5».  ^  ^  2 
(0  —  a)Uy  —  (1  -  p)s 


1  —  i?043  —  Cl(/l) 


(A.18) 


Rim  -  Ro43  =  fc  +  o(fc2)  ^  0(^2, , 

(b  —  a)uy  —  (1  —  p)s 


Then  because  is  Lipschitz  continuous 

(A.19)  l-«'(i2o43)  =  (l -«'(!)) +  0(/i) 
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and 

(A. 20)  ^(jRyso)  —  'J'(.Ro43)  =  0{h^), 

which  together  with  the  fact  that  the  2D  scheme  is  second  order  accurate  proves  the  lemma.  □ 


Appendix  B.  Fluctuation-splitting  scheme  for  the  inhomogeneous  case 

Here  we  address  the  problem  of  maintaining  linearity  preserving  property  the  NNL  scheme  on 
unstructured  triangular  meshes  in  the  case  of  non-zero  forcing  term.  Consider  again  the  inhomo¬ 
geneous  advection  equation  (4.1).  The  N  scheme  given  by  (4. 16), (4. 17)  for  the  homogeneous  case 
should  be  now  modified.  One  target  formula  (Type  I  triangle,  Fig. 7)  now  is 

(B.l)  =  Szul  +  At[ki{ui  -  Us)  +  k2{u2  -  us)  -  cxt] 

and  in  two  target  case  (Type  II  triangle,  Fig.7) 

=  Siu'l;  +  At[ki{ui  -  Us)  -  fiaj] 

^  ’  S2U'^+^  =  S2U^  +  At[k2{u2  -  Us)  -  (1  -  P)aTl 

where 


^-kr  +  k2 

and  (Tt  defined  by  (4.4).  The  following  distribution  formulae  will  give  the  inhomogeneous  NNL 
scheme  in  the  two  target  case 

5iu^+i  =  5iUi  -h  At{[  ki  (ui  -  Us)  -  fiar] 

.X  +'^{Rl32)[k2iu2  -Us)  -{I-  P)(Tt]} 

^  S 2U2'^^  =  Ssu'^  +  At {[  k2{u2  -  Us)  -  {1  -  l3)aT] 

-'i{Rl32)[k2{U2  -  Us)  -  (1  -  P)(rT]}, 


where 


Identity 


ki{ui  -  Us)  -  l3aT 
k2{u2  -  Us)  -  (1  -  /3)c7-r 


(B.6)  ^{Ri32)[k2{u2  -  Us)  -  (1  -  I3)(tt]  =  -  us)  -  Par] 

til32 

Positivity  of  this  scheme  can  be  demonstrated  in  the  same  way  as  for  the  homogeneous  equation. 
The  following  result  concerning  the  LP  property  of  the  constructed  scheme  can  be  formulated 

Theorem  B.l.  The  NNL  scheme  (B.1),(B.2)  is  linearity  preserving  in  the  inhomogeneous  case  if 

(B.7)  |'F(i?)-l|<C'|fi-li 

for  some  positive  constant  C. 

Proof.  Define 

Sui  =  Siu^'^^  -  Siu'^ 

=  At{[ki{ui  -  Us)  -  Pc^t]  +  '9{Rl32)[k2{u2  -  Us)  -  (1  -  P)crT]} 

=  P  +  ^{R)Q 

(B.8) 

5u2  =  SsU^'^^  -  S2U2 

=  At{[k2{u2  -  Us)  -  (1  -  P)crT]  -  {Rl32)[k2{U2  -  Us)  -  (1  -  P)o-t]} 

=  Q-^{R)Q. 

The  rest  of  the  proof  is  the  same  as  of  Theorem  5.3.  □ 


To  demonstrate  the  connection  between  the  inhomogeneous  NNL  scheme  and  previously  con¬ 
structed  S*  scheme  consider  the  grid  illustrated  in  Fig. 9.  The  fluctuations  from  the  triangles  A,  B,  C 
can  be  given  by  the  following  formulae 

[1  -  ^(i?043)][(6-  a)(«3  -  «4)  -  (1  -  ^)hs^] 

+a(uo  -  Uj)  -  +  fl*(fio43)[(^-  g)('^3  -  ^4)  -  (1  -  f^)hs^]} 

(B.9)  = -fi  q(^o  -  •^4)  +  (h  -  a){uo  -  U5)  -  h^} 

=  -fi  -  Uz)  -  +  ^(i?75o)[(6  -  a)(rio  -  ^s)  -  (1  -  P)hs'^] 

+  [1  -  ’I>(i?75o)][(^-  «)(^0  -  ^s)  -  (1  -  li)hs‘^]), 

where  the  underlined  terms  contribute  to  the  residual  at  the  node  0.  Thus  the  update  formula  is 
=Uq  +  ^{  b{u5-uo)  +  a{u4  -  U5) 

-^fl'(i?75o)[(i  -  a){n5  -  uo)  +  h{l  -  f3)sc] 

-t-^fl'(/?043)[(h-  a){u4  -  Us)  +  h{l  -  (3)sa] 

(B.IO)  +~[sb  +  (1  +  ^)<5.4  +  (1  ~ /?)'5c]}- 

Recalling  the  definition  of  sa,sb,sc  it  is  easy  to  see  that  (B.IO)  is  identical  to  (A. 10). 

ICASE,  Mail  Stop  132C,  NASA  Langley  Research  Center,  Hampton,  VA  23681 
E-mail  address:  sidilkov00icase.edu 

Department  of  Aerospace  Engineering,  University  of  Michigan,  Ann  Arbor,  MI  48109. 

E-mail  address:  philroe00engin.uittich.edu 


26 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  includingthe  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson 
Davis  Highway,  Suite  1204,  Arlington.  VA  22202-4302.  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188).  Washington,  DC  20503. 

1.  AGENCY  USE  ONLY (Uave  blank) 

2.  REPORT  DATE 

February  1995 

3.  REPORT  TYPE  AND  DATES  COVERED 

Contractor  Report 

4.  TITLE  AND  SUBTITLE 

UNIFICATION  OF  SOME  ADVECTION  SCHEMES  IN  TWO 
DIMENSIONS 

5.  FUNDING  NUMBERS 

C  NASl-19480 

WU  505-90-52-01 

6.  AUTHOR(S) 

D.  Sidilkover 

P.  L.  Roe 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Institute  for  Computer  Applications  in  Science 
and  Engineering 

Mail  Stop  132C,  NASA  Langley  Research  Center 

Hampton,  VA  23681-0001 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

ICASE  Report  No.  95-10 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

N  ational  Aeronautics  and  Space  Administration 

Langley  Research  Center 

Hampton,  VA  23681-0001 

10.  SPONSORING/MONITORING 

AGENCY  REPORT  NUMBER 

NASA  CR-195044 

ICASE  Report  No.  95-10 

n.  SUPPLEMENTARY  NOTES 

Langley  Technical  Monitor:  Dennis  M.  Bushnell 

Final  Report 

Submitted  to  Math.  Comp. 

12a.  DISTRIBUTION/AVAILABILITY  STATEMENT 

U  nclassified-Unlimited 

Subject  Category  64 

12b.  DISTRIBUTION  CODE 

13.  ABSTRACT  (Maximum  200  wards) 

In  this  paper  a  relationship  between  two  approaches  towards  construction  of  genuinely  two-dimensional  upwind 
advection  schemes  is  established.  One  of  these  approaches  is  of  the  control  volume  type  applicable  on  structured 
cartesian  meshes.  It  resulted  (see  [14],  [15])  in  the  compact  high  resolution  schemes  capable  of  maintaining  second 
order  accuracy  in  both  homogeneous  and  inhomogeneous  cases.  Another  one  is  the  fluctuation  splitting  approach 
(see  [11],  [3],  [12],  [17]),  which  is  well  suited  for  triangular  (and  possibly)  unstructured  meshes.  Understanding  the 
relationship  between  these  two  approaches  allows  us  to  formulate  here  a  new  fluctuation  splitting  high  resolution 
(i.e.  possible  use  of  artificial  compression,  while  maintaining  positivity  property)  scheme.  This  scheme  is  shown  to 
be  linearity  preserving  in  inhomogeneous  as  well  as  homogeneous  cases. 

14.  SUBJECT  TERMS 

multidimensional  advection;  non-zero  forcing  term;  structural  and  unstructural  meshes 

15.  NUMBER  OF  PAGES 

28 

16.  PRICE  CODE 

A03 

17.  SECURITY  CLASSIFICATION 

OF  REPORT 

Unclassified 

18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

Unclassified 

19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

20.  LIMITATION 

OF  ABSTRACT 

NSN  7540-01-280-5500  Standard  Form  298(Rev.  2-89) 

Prescribed  by  ANSI  Std.  Z39-18 


298-102 


